/// <license>
/// This is a port of the SciMark2a Java Benchmark to C# by
/// Chris Re (cmr28@cornell.edu) and Werner Vogels (vogels@cs.cornell.edu)
/// 
/// For details on the original authors see http://math.nist.gov/scimark2
/// 
/// This software is likely to burn your processor, bitflip your memory chips
/// anihilate your screen and corrupt all your disks, so you it at your
/// own risk.
/// </license>

using System;

namespace SpeakerRecognition
{
	
	/// <summary>Computes FFT's of complex, double precision data where n is an integer power of 2.
	/// This appears to be slower than the Radix2 method,
	/// but the code is smaller and simpler, and it requires no extra storage.
	/// </P>
	/// </summary>
	/// 
	/// <author> 
	/// Bruce R. Miller bruce.miller@nist.gov,
	/// Derived from GSL (Gnu Scientific Library), 
	/// GSL's FFT Code by Brian Gough bjg@vvv.lanl.gov
	/// </author>
	
	public class FFT
	{
        // compute the FFT of x[], assuming its length is a power of 2
        public static Complex[] fft(Complex[] x)
        {
            int N = x.Length;
            Complex[] y = new Complex[N];

            // base case
            if (N == 1)
            {
                y[0] = x[0];
                return y;
            }

            // radix 2 Cooley-Tukey FFT
            if (N % 2 != 0) throw new Exception("N is not a power of 2");
            Complex[] even = new Complex[N / 2];
            Complex[] odd = new Complex[N / 2];
            for (int k = 0; k < N / 2; k++) even[k] = x[2 * k];
            for (int k = 0; k < N / 2; k++) odd[k] = x[2 * k + 1];

            Complex[] q = fft(even);
            Complex[] r = fft(odd);

            for (int k = 0; k < N / 2; k++)
            {
                double kth = -2 * k * Math.PI / N;
                Complex wk = new Complex(Math.Cos(kth), Math.Sin(kth));
                y[k] = q[k].plus(wk.times(r[k]));
                y[k + N / 2] = q[k].minus(wk.times(r[k]));
            }
            return y;
        }

        // compute the inverse FFT of x[], assuming its length is a power of 2
        public static Complex[] ifft(Complex[] x)
        {
            int N = x.Length;
            Complex[] y = new Complex[N];

            // take conjugate
            for (int i = 0; i < N; i++)
            {
                y[i] = x[i].conjugate();
            }

            // compute forward FFT
            y = fft(y);

            // take conjugate again
            for (int i = 0; i < N; i++)
            {
                y[i] = y[i].conjugate();
            }

            // divide by N
            for (int i = 0; i < N; i++)
            {
                y[i] = y[i].times(1.0 / N);
            }

            return y;

        }

        public static Complex[] convertToComplex(double[] arr)
        {
            Complex[] x = new Complex[arr.Length];
            // original data
            for (int i = 0; i < arr.Length; i++)
                x[i] = new Complex(arr[i], 0);

            return x;
        }
        /// <summary>
        /// findMin to calculate the magnitudes of a complex number using sqrt of the sum of the powers
        /// if the realOnly flag is true then the abs real part 
        /// </summary>
        /// <param name="input"></param>
        /// <param name="realOnly">realOnly flag for abs real part </param>
        /// <returns></returns>
        public static double[] getMagnitude(Complex[] input, bool realOnly)
        {
            double[] magnitudeData = new double[input.Length];
            if (realOnly)
            {
                for (int i = 0; i < input.Length; i++)
                {
                    magnitudeData[i] = Math.Sqrt(input[i].real() * input[i].real());
                }
            }
            else
            {
                for (int i = 0; i < input.Length; i++)
                {
                    magnitudeData[i] = input[i].mod();
                }
            }
            return magnitudeData;
        }

	}
}